Análisis espacio-temporal datos histórica artesanal de la pesquería de Bacalao

Analisis exploratorio 1986-2021

Author

Mauricio Mardones I

Published

July 13, 2023

ANTECEDENTES

El presente reporte tiene un codigo autocontenido con los Analísis Exploratorios de Datos (AED) de la pesquería de bacalao de profundidad extraído por la flota pesquera artesanal (espinel). Los analisis aca descritos estan compuestos de dos bases. La primera es de los registros oficiales de captura recogidos por el SERNAPESCA. En segundo lugar, se presentan el AED de los datos del monitoreo d la pesqería llevado a cabo por el IFOP a través del Programa de Seguimiento de Pesquerías Demersales del Departamento de Evaluación de Pesquerías.

El principal objetivo es identificar vacios y fortalezas de los datos en terminos espacio-temporales, y con ello tomar decisiones para el paso posterior de modelación de la dinámica poblacional del recurso para realizar asesoría en terminos de estatus y recomendación de una Captura Biologicamente Aceptable (CBA).

ANÁLISIS EXPLORATORIO DE DATOS (AED)

Cargo librerías necesarias para el análisis exploratorio de los datos de las distintas bases de bitácora, tallas y biológico.

library(here)
#analisis
library(ggsignif)
library(ggrepel)
library(ggpubr)
library(inlmisc)
library(nortest) #para testear distribucion
library(skimr) #provides a frictionless approach to summary statistics w
library(easystats) # multiples unciones analiticas
library(lme4)
library(skimr)
library(readxl)
# vizualizacion
library(ggridges)
library(sf)
library(GGally)
library(tidyverse, quietly = TRUE)
library(knitr, quietly = TRUE)
library(kableExtra)
library(raster)
library(egg)
library(car) #Variance inflation Factor
library(ggthemes)
library(sjPlot)
# A function for dotplots
multi_dotplot <- function(filename, Xvar, Yvar){
  filename %>%
    ggplot(aes(x = {{Xvar}})) +
    geom_point(aes(y = {{Yvar}}),
               alpha=0.4) +
    theme_bw() +
    coord_flip() +
    labs(x = "Order of Data")}

Identifico los directorio de trabajo y leo las tres bases de datos; a saber:

  • Datos FIP (Rendimiento) 96-32
  • Bitácora
  • Biologico
  • Tallas
  • Desembarques (1985-2022)
bit <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BITACORAS ESPINEL bacalao_corr.xlsx", sheet = "BITACORAS_ESPINEL_bacalao",  guess_max = 100000)
long <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/LONGITUD_ESPINEL_bacalao.xlsx")
bio <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BIOLOGICO ESPINEL bacalao.xlsx", 
    sheet = "BIOLOGICO_ESPINEL_bacalao")
landing <- read_excel("Artesanal_Historica/DESEMBARQUE HISTÓRICO.xlsx", 
    skip = 2)
bit8696<- read_csv2("bit_art_86_20.csv")

Identificamos los registros asociados recurso objetivo

  • Codigo del Recurso: 37
  • Nombre común: Bacalao de profundidad
  • Nombre científico: Dissostichus eleginoides

Idetifico los registros por recurso y por base para luego filtrar.

dim(bit)
   [1] 14208    77
table(bit$COD_ESPECIE)
   
       2     4     5     6    14    15    16    17    18    20    22    23    24 
       7     1    14    25   147     6    89    16   698     1    25     3     7 
      25    27    33    35    37    42    43    49    53    56    59    74    75 
       3    11     2     2 11074     3     9    21     1     5    35     1    28 
      81    88    99   103   104   105   106   108   119   127   129   144   194 
      19     1   159    21    24     2     7    44     1     1     1     1    91 
     199   212   226   239   243   323   324   325   334   844   849   920   984 
     178    18     1    22     1     2    50     6   830     3   138     1     1 
     989   999  1309 
      36   249    58
dim(long)
   [1] 15831    24
table(long$COD_ESPECIE)
   
      37 
   15831
dim(bio)
   [1] 127409     36
table(bio$COD_ESPECIE)
   
       37 
   127409
# las bases de biologicos y longitud solo teienen  registros de bacalao.

Comenzar a trabajar bases por separado

DESEMBARQUES

Tabla con los desembarques oficiales (Sernapesca, 2022)


kbl(landing, booktabs = T,format = "html",
    caption = "Desembarque Bacalao Artesanal por Región") %>%
    kable_styling(latex_options = c("striped",
                                  "condensed","scale_down"),
                full_width = FALSE) 
Desembarque Bacalao Artesanal por Región
...1 AyP TPCA ANTOF ATCMA COQ VALPO LGBO MAULE ÑUBLE BBIO ARAUC RIOS LAGOS AYSEN MAG ...17
1985 NA 8 102 146 162 1674 NA 349 NA 1527 NA NA 63 28 NA 4059
1986 NA 842 1531 213 329 1242 NA 88 NA 1546 NA NA 311 6 NA 6108
1987 NA 182 335 171 51 671 NA NA NA 1264 NA NA 689 21 NA 3384
1988 NA 131 174 251 42 670 NA 46 NA 1597 NA NA 877 8 NA 3796
1989 NA 217 175 177 109 756 NA 114 NA 1859 10 NA 1462 8 NA 4887
1990 NA 335 161 126 100 444 NA 5 NA 2620 NA NA 1789 36 NA 5616
1991 NA 280 186 142 223 519 NA NA NA 1901 NA NA 680 NA NA 3931
1992 NA 30 45 77 53 477 NA 131 NA 2108 1 NA 741 1 NA 3664
1993 NA 37 4 66 81 1676 NA 295 NA 1111 NA NA 807 NA 45 4122
1994 NA 24 139 213 68 671 NA 526 NA 1920 NA NA 1374 41 411 5387
1995 NA 252 259 286 65 206 NA 190 NA 1347 NA NA 1680 31 266 4582
1996 NA 476 220 121 177 272 NA 226 NA 1390 2 NA 1989 22 83 4978
1997 NA 36 82 47 57 98 NA 116 NA 1077 1 NA 1803 22 83 3422
1998 NA 9 67 169 61 155 NA 399 NA 1284 NA NA 2038 NA 11 4193
1999 NA 42 72 168 87 484 NA 406 NA 1369 NA NA 2910 NA 270 5808
2000 NA 285 98 153 105 200 NA 292 NA 800 NA NA 3515 NA 345 5793
2001 NA 198 85 127 113 191 NA 136 NA 1049 NA NA 2045 NA NA 3944
2002 NA 154 44 97 57 188 NA 186 NA 732 NA NA 3097 NA 10 4565
2003 NA 58 25 51 27 338 NA 123 NA 752 NA NA 3189 NA 179 4742
2004 NA 76 32 63 23 308 NA 132 NA 412 NA NA 2303 NA 70 3419
2005 NA 111 14 47 26 84 NA 175 NA 475 NA NA 2021 NA 325 3278
2006 NA 87 3 50 50 92 NA 159 NA 243 NA NA 1367 NA 40 2091
2007 19 43 7 22 25 43 NA 84 NA 189 NA 501 1157 NA NA 2090
2008 4 48 21 9 3 32 NA 70 NA 219 NA 349 803 NA NA 1558
2009 13 85 24 20 16 79 NA 112 NA 239 NA 522 571 NA NA 1681
2010 13 43 13 58 11 37 NA 63 NA 109 NA 465 655 NA NA 1467
2011 19 82 20 46 6 49 NA 90 NA 276 NA 613 987 NA NA 2188
2012 25 68 11 82 14 26 NA 105 NA 317 NA 290 1126 NA NA 2064
2013 26 58 12 68 23 52 NA 214 NA 186 NA 277 622 NA 20 1558
2014 13 48 - 59 34 64 NA 171 NA 129 NA 252 453 NA 57 1280
2015 6 52 13 53 38 93 NA 119 NA 261 NA 298 638 NA 38 1609
2016 17 51 20 59 16 70 NA 71 NA 335 NA 327 755 NA NA 1721
2017 20 92 33 89 33 117 NA 116 NA 316 NA 312 790 NA 28 1946
2018 29 78 5 74 49 78 NA 113 NA 251 NA 218 620 NA 25 1540
2019 26 55 7 63 35 42 NA 102 NA 248 NA 203 696 NA 94 1571
2020 5 18 7 NA 1 27 NA 78 NA 145 NA 229 261 NA 32 803
2021 2 83 10 138 17 52 NA 206 NA 304 NA 500 517 NA 21 1850
2022 26 100 - 148 58 75 - 143 - 318 - 332 754 - 29 1983

Ploteo los desembarques por region y por año

landing <- as.data.frame(lapply(landing, as.double))
str(landing)
   'data.frame':    38 obs. of  17 variables:
    $ ...1 : num  1985 1986 1987 1988 1989 ...
    $ AyP  : num  NA NA NA NA NA NA NA NA NA NA ...
    $ TPCA : num  8 842 182 131 217 335 280 30 37 24 ...
    $ ANTOF: num  102 1531 335 174 175 ...
    $ ATCMA: num  146 213 171 251 177 126 142 77 66 213 ...
    $ COQ  : num  162 329 51 42 109 100 223 53 81 68 ...
    $ VALPO: num  1674 1242 671 670 756 ...
    $ LGBO : num  NA NA NA NA NA NA NA NA NA NA ...
    $ MAULE: num  349 88 NA 46 114 5 NA 131 295 526 ...
    $ ÑUBLE: num  NA NA NA NA NA NA NA NA NA NA ...
    $ BBIO : num  1527 1546 1264 1597 1859 ...
    $ ARAUC: num  NA NA NA NA 10 NA NA 1 NA NA ...
    $ RIOS : num  NA NA NA NA NA NA NA NA NA NA ...
    $ LAGOS: num  63 311 689 877 1462 ...
    $ AYSEN: num  28 6 21 8 8 36 NA 1 NA 41 ...
    $ MAG  : num  NA NA NA NA NA NA NA NA 45 411 ...
    $ ...17: num  4059 6108 3384 3796 4887 ...
summary(landing)
         ...1           AyP             TPCA            ANTOF        
    Min.   :1985   Min.   : 2.00   Min.   :  8.00   Min.   :   3.00  
    1st Qu.:1994   1st Qu.:11.25   1st Qu.: 44.25   1st Qu.:  12.75  
    Median :2004   Median :18.00   Median : 77.00   Median :  32.50  
    Mean   :2004   Mean   :16.44   Mean   :128.26   Mean   : 112.67  
    3rd Qu.:2013   3rd Qu.:25.25   3rd Qu.:148.25   3rd Qu.: 111.25  
    Max.   :2022   Max.   :29.00   Max.   :842.00   Max.   :1531.00  
                   NA's   :22                       NA's   :2        
        ATCMA            COQ             VALPO             LGBO    
    Min.   :  9.0   Min.   :  1.00   Min.   :  26.0   Min.   : NA  
    1st Qu.: 58.0   1st Qu.: 23.50   1st Qu.:  65.5   1st Qu.: NA  
    Median : 82.0   Median : 49.50   Median : 136.0   Median : NA  
    Mean   :106.7   Mean   : 64.34   Mean   : 325.1   Mean   :NaN  
    3rd Qu.:148.0   3rd Qu.: 77.75   3rd Qu.: 468.8   3rd Qu.: NA  
    Max.   :286.0   Max.   :329.00   Max.   :1676.0   Max.   : NA  
    NA's   :1                                         NA's   :38   
        MAULE           ÑUBLE          BBIO            ARAUC           RIOS      
    Min.   :  5.0   Min.   : NA   Min.   : 109.0   Min.   : 1.0   Min.   :203.0  
    1st Qu.: 99.0   1st Qu.: NA   1st Qu.: 253.5   1st Qu.: 1.0   1st Qu.:270.8  
    Median :127.0   Median : NA   Median : 603.5   Median : 1.5   Median :319.5  
    Mean   :165.3   Mean   :NaN   Mean   : 848.0   Mean   : 3.5   Mean   :355.5  
    3rd Qu.:194.0   3rd Qu.: NA   3rd Qu.:1363.5   3rd Qu.: 4.0   3rd Qu.:473.8  
    Max.   :526.0   Max.   : NA   Max.   :2620.0   Max.   :10.0   Max.   :613.0  
    NA's   :2       NA's   :38                     NA's   :34     NA's   :22     
        LAGOS            AYSEN            MAG             ...17     
    Min.   :  63.0   Min.   : 1.00   Min.   : 10.00   Min.   : 803  
    1st Qu.: 661.2   1st Qu.: 8.00   1st Qu.: 28.25   1st Qu.:1753  
    Median : 842.0   Median :22.00   Median : 51.00   Median :3402  
    Mean   :1267.2   Mean   :20.36   Mean   :112.82   Mean   :3228  
    3rd Qu.:1799.5   3rd Qu.:29.50   3rd Qu.:157.75   3rd Qu.:4472  
    Max.   :3515.0   Max.   :41.00   Max.   :411.00   Max.   :6108  
                     NA's   :27      NA's   :16
landing2 <- landing %>% 
  pivot_longer(cols=c("AyP"  , "TPCA"  ,"ANTOF", "ATCMA", "COQ",
               "VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
               "RIOS" , "LAGOS", "AYSEN", "MAG"), 
               names_to = "REGION", 
               values_to = "CAPTURA") %>% 
  rename(AÑO="...1") %>% 
  dplyr::select(-2)

Genero un gráfico de barras


landing2$REGION <- factor(landing2$REGION , 
                          levels = c("AyP"  , "TPCA"  ,"ANTOF", "ATCMA", "COQ",
               "VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
               "RIOS" , "LAGOS", "AYSEN", "MAG"))
desem <- ggplot(landing2 %>% 
                   drop_na(REGION),aes(AÑO, CAPTURA, fill=REGION)) +
  geom_bar(stat="identity")+
  scale_fill_viridis_d(option="G")+
  theme_minimal()+
  scale_x_continuous(breaks = seq(from = 1985, to = 2022, by = 4))+
  scale_y_continuous(breaks = seq(from = 0, to = 4000, by = 1000))+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank(),
        legend.position = "none")+
  facet_wrap(~REGION, ncol=3)+
  labs(y="Captura (t)",
       x="")
desem

Los principales desembarques estan asociados a las regiones de Antifagasta, Valparaiso y BioBIo pero solo durante los primeros años (1985- mediados del 2000). Luego de esto, todas las regiones vieron disminuidos los registros de extracción del recurso.

BITÁCORA

Esta base tiene como principal objetivo obtener un indicador de esfuerzo para calcular un indice de abundancia relativo como la CPUE.

Primero identifico la estructura de la base y filtro el rrecurso bacalao

# filtro bitacoras dejando solo bacalao
bitb <- bit %>% 
  filter(COD_ESPECIE==37)%>% 
  mutate(dfp = as.double(dfp), 
         PESO_corr = as.double(PESO_corr),
         prof_med=as.double(prof_med))
dim(bitb)
   [1] 11074    77

se eliminan 1208-11074 = 3134 registros

glimpse(bitb)

ahora las estadísticas descriptivas de las variables de intéres. En este caso dfp, PESO_corr, prof_med

Miramos los oultiers de los datos

# OUTLIERS

#Order data
bitb <- bitb %>%
  mutate(order = seq(1:nrow(bitb)))

#Select continuous variables to plot
p1 <- multi_dotplot(bitb, order, dfp)
p2 <- multi_dotplot(bitb, order, PESO_corr)
p3 <- multi_dotplot(bitb, order, prof_med)
p4 <- multi_dotplot(bitb, order, N_TRIPULANTES)
p5 <- multi_dotplot(bitb, order, HORAS_REPOSO)
p6 <- multi_dotplot(bitb, order, NUMERO_DE_ANZUELOS)

#Plot as a grid
ggarrange(p1, p2, p3, p4,p5, p6,  nrow = 3)

Remuevo outliers

Los datos de profundidades andan bien, Filtro los datos de la variable dfp entre 0 y 30 días. y repito la inspección

bitb2 <- bitb  %>% 
  filter(dfp>0,
         dfp<80,
         PESO_corr<50000,
         N_TRIPULANTES<30,
         NUMERO_DE_ANZUELOS<40000)

Miro nuenvamente los datos filtrados

# OUTLIERS

#Order data
bitb2 <- bitb2 %>%
  mutate(order = seq(1:nrow(bitb2)))

#Select continuous variables to plot
p1a <- multi_dotplot(bitb2, order, dfp)
p2a <- multi_dotplot(bitb2, order, PESO_corr)
p3a <- multi_dotplot(bitb2, order, prof_med)
p4a <- multi_dotplot(bitb2, order, N_TRIPULANTES)
p5a <- multi_dotplot(bitb2, order, HORAS_REPOSO)
p6a <- multi_dotplot(bitb2, order, NUMERO_DE_ANZUELOS)

#Plot as a grid
ggarrange(p1a, p2a, p3a, p4a, p5a, p6a,  nrow = 3)

Identifico las dimensionesd e la nueva base y la comparo con la anterior

dim(bitb2)
   [1] 4745   78
dim(bitb)
   [1] 11074    78
table(bitb2$año)
   
   2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 
      2    3   59  396  441  364  754  750  468  497  561  450

Los registros se minimiza en ciertos años lo cual genera problemas en la fiabilidad de la estimación.

ahora calculamos el rendimiento nominal

CPUE Nominal

Ahora estimo la CPUE, la cual es el rendimoento entre el esfuerzo y la captura. La medida de esfuerzo será dfp. También calculo un rendimiento con la variable de esfuerzo HORAS_REPOSO pero solo como antecedente dado que no tengo muchos datos.

bitb2 <- bitb2 %>% 
  mutate(CPUE = PESO_corr/dfp) %>% 
  mutate(CPUE2 =PESO_corr/HORAS_REPOSO)
bitb2$prof_med_cat <- cut(bitb2$prof_med, 
                          breaks = c(500, 1000, 1500,
                                     2000, 3000))
# Frequency polygon plot for catch
his <- ggplot(bitb2, aes(dfp)) +
  geom_freqpoly(bins = 5) +
  labs(x = "dfp", y = "Frequency") +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", 
                                    fill=NA, size = 1))
# Patterns in the variance? (any lack of homogeneity)
point <- ggplot(bitb2 %>% 
                  filter(PESO_corr>1000), aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 5, alpha = 0.6) +
  stat_smooth(method = "lm")+
  theme(panel.background = element_blank()) +
  theme(panel.border = element_rect(fill = NA, size = 1)) +
  theme(strip.background = element_rect(fill = "white", 
                   color = "white", size = 1)) +
  theme(text = element_text(size=13)) +
  xlab("dfp") + ylab("Peso Corregido")


#Plot as a grid
ggarrange(his, point,  nrow = 1)

Tendencia del esfuerzo dfp a través de los años y por región

effor <- ggplot(bitb2 %>%
                  group_by(año, REGION_PUERTO_RECALADA) %>% 
                  summarise(DFPM=mean(dfp)),
                aes(año, DFPM))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
  geom_smooth(method = 'lm', 
              colour = 'blue', 
              size = 1.5)+
  facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="Effort (dfp)",
       x="")
effor

Tendencia del esfuerzo NUMERO_DE_ANZUELOS a través de los años y por región

efforanz <- ggplot(bitb2 %>%
                  group_by(año, REGION_PUERTO_RECALADA) %>% 
                  summarise(NANZ=mean(NUMERO_DE_ANZUELOS)),
                aes(año, NANZ))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
  geom_smooth(method = 'lm', 
              colour = 'green', 
              size = 1.5)+
  facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="Effort (dfp)",
       x="")
efforanz

Calculo los ceros de dfp

#CALCULATE NUMBER OF ZEROS

# What is the percentage of zeros i the response variable

round(sum(bitb2$dfp == 0) * 100 / nrow(bitb2),0)
   [1] 0
#0

Interacciones

# Interactions

# Year x season
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("Points") + ylab("Catch") +
  facet_grid(año~trim)

# No

ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~N_TRIPULANTES)

# Perhaps

ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(trim~N_TRIPULANTES)

# CPUE slope varies between habitats - interaction

Profundidad por region

profu <- ggplot(bitb2, aes(REGION_PUERTO_RECALADA, desc(prof_med), 
                            group=REGION_PUERTO_RECALADA))+
  geom_boxplot(fill=2, alpha=.5)+
  geom_jitter(size=0.4, alpha=0.2,
              width = .25)+
  #facet_wrap(~año)+
  theme_bw()+
  scale_x_continuous(breaks = seq(from = 1, to = 15, by = 1))+
  labs(x="Región",
       y="Profundidad (mts)")
profu

ggplot(bitb2, aes(x = nombre_puerto_recalada , y = dfp, fill = nombre_puerto_recalada )) +
  geom_violindot(fill_dots = "black") +
  theme_modern() +
  scale_fill_material_d()


ggplot(bitb2, aes(x = as.factor(REGION_PUERTO_RECALADA ), 
                  y = bitb2$NUMERO_DE_ANZUELOS, 
                  fill = as.factor(REGION_PUERTO_RECALADA))) +
  geom_violin() +
  theme_modern() +
  scale_fill_material_d(palette="ice",
                        name="REGION")

ggplot(bitb2 %>% 
         drop_na(prof_med_cat), 
       aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', 
              colour = 'red', 
              se = FALSE) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(prof_med_cat~año)

ploteo los datos de CPUE totales.

cpuenom <- ggplot(bitb2 %>% 
  group_by(zona, año) %>%
  summarise(CPUEM=mean(CPUE)) %>% 
  filter(CPUEM<600,
         año>1996,
         zona!=4),
  aes(año, CPUEM, color=zona, group=zona))+
    geom_point()+
  geom_smooth()+
  scale_color_viridis_b(option="B",
                        labels = seq(1, 3, by = 1))+
  scale_x_continuous(breaks = seq(from = 2010, to = 2022, by = 1))+
  theme_few()+
  labs(x="",
       y="CPUE")
cpuenom

guardo los datos bitacora CPUE nominal 1996-2022 del objeto cpuenom

write_csv(cpuenom, "CPUE_MON.csv")

CPUE Estandarizada

Identificar los proncipales factores para modelar la variable.

Luego de una reunión con Seguimiento (Patricio Galvez) se indicaron algunos aspectos que deben ser considerados en la estandariación.

  • El poder de pesca (tipo de embarcación) es mayor en la zona centro sur de Chile. Identificar equilibrio de la base bitb2 respecto a este dato.

  • Usar zonas como factor en desmedro de región.

Identifico una estadistica descriptiva general de la base de estandarización

Correlaciones

Análisis de utilidad para identificar a traves de un metodo de corrrelación de pearson, la correlación en tre las variables que serán utilizadas en la estandariacion de la CPUE.

Primero identifico las variables a correlacionar

names(bitb2)
    [1] "ID"                           "embarque"                    
    [3] "COD_BARCO"                    "FECHA_HORA_RECALADA"         
    [5] "FECHA_HORA_ZARPE"             "COD_PESQUERIA"               
    [7] "COD_PESQUERIA_ant"            "N_TRIPULANTES"               
    [9] "N_CALADAS"                    "PUERTO_ZARPE"                
   [11] "PUERTO_RECALADA"              "nombre_puerto_recalada"      
   [13] "NRO_FORMULARIO"               "REGION_PUERTO_RECALADA"      
   [15] "ESPECIE_OBJETIVO"             "ECOSONDA"                    
   [17] "GPS"                          "POTENCIA_MOTOR_EQ"           
   [19] "VIRADOR"                      "GASTOS_VIAJE"                
   [21] "NUMERO_LANCE_EX"              "FECHA_LANCE"                 
   [23] "año"                          "mes"                         
   [25] "trim"                         "dfp"                         
   [27] "ID_CUADRICULA"                "ID_PROCEDENCIA"              
   [29] "PESO_TOTAL_CAPTURA"           "LATITUD"                     
   [31] "LONGITUD"                     "lat_ctgr"                    
   [33] "lon_ctgr"                     "zona"                        
   [35] "ESPECIE_OBJETIVO_LANCE"       "CLASE_LANCE"                 
   [37] "COD_ESPECIE"                  "PESO"                        
   [39] "PESO_corr"                    "CAPTURA_RETENIDA"            
   [41] "VOLUMEN_DESCARTE"             "LUGAR_DESCARTE"              
   [43] "TIPO_DESCARTE"                "PRECIO_UNITARIO"             
   [45] "DESTINO_CAPTURA"              "NRO_ANZUELOS"                
   [47] "REGION_PROCEDENCIA"           "FECHA_HORA_FIN_CALADO"       
   [49] "FECHA_HORA_INI_CALADO"        "LATITUD_INICIO_CALADO"       
   [51] "LONGITUD_INICIO_CALADO"       "LATITUD_FIN_CALADO"          
   [53] "LONGITUD_FIN_CALADO"          "FECHA_HORA_INI_VIRADO"       
   [55] "FECHA_HORA_FIN_VIRADO"        "LATITUD_INICIAL_VIRADO"      
   [57] "LONGITUD_INICIAL_VIRADO"      "LATITUD_FIN_VIRADO"          
   [59] "LONGITUD_FIN_VIRADO"          "PROFUNDIDAD_MINIMA_LM"       
   [61] "PROFUNDIDAD_MAXIMA_LM"        "PROFUNDIDAD_FONDO_VIRADO_INI"
   [63] "PROFUNDIDAD_FONDO_VIRADO_FIN" "prof_med"                    
   [65] "SEPARACION_ANZUELO"           "MARCA_ANZUELO"               
   [67] "NUMERO_DE_ANZUELOS"           "PERDIDA_ANZUELOS"            
   [69] "TIPO_CARNADA"                 "TAMANIO_ANZUELO"             
   [71] "INTENSIDAD_VIENTO"            "LONGITUD_LINEA_MADRE"        
   [73] "HORAS_REPOSO"                 "TIPO_ESPINEL"                
   [75] "N_PANOS"                      "ANZ_POR_PANO"                
   [77] "filtro"                       "order"                       
   [79] "CPUE"                         "CPUE2"                       
   [81] "prof_med_cat"
bitb3 <- bitb2 %>%
  dplyr::select(8, 18, 23,24, 25, 26, 32, 33, 34, 39, 69, 73, 75, 78, 79, 80, 81)
results <- correlation(bitb3)

results %>%
  summary(redundant=TRUE) %>% 
  plot(results, show_data = "points")+
  theme_bw()

skimr::skim(bitb3)
Data summary
Name bitb3
Number of rows 4745
Number of columns 17
_______________________
Column type frequency:
character 1
factor 1
numeric 15
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
TIPO_CARNADA 122 0.97 1 5 0 8 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
prof_med_cat 939 0.8 FALSE 4 (1e: 1867, (50: 961, (1.: 853, (2e: 125

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
N_TRIPULANTES 0 1.00 7.23 1.62 2.00 6.00 7.00 9.00 10.00 ▁▂▆▇▆
POTENCIA_MOTOR_EQ 167 0.96 405.74 399.03 0.00 320.00 390.00 420.00 3903.00 ▇▁▁▁▁
año 0 1.00 2018.05 2.47 2011.00 2016.00 2018.00 2020.00 2022.00 ▁▅▆▆▇
mes 0 1.00 7.29 3.79 1.00 3.00 9.00 11.00 12.00 ▅▂▁▂▇
trim 0 1.00 2.71 1.32 1.00 1.00 3.00 4.00 4.00 ▆▂▁▂▇
dfp 0 1.00 16.73 6.62 0.67 13.02 16.79 19.79 50.67 ▂▇▂▁▁
lat_ctgr 51 0.99 -35.49 6.89 -55.17 -40.99 -35.25 -33.17 -18.33 ▁▅▇▂▂
lon_ctgr 51 0.99 -73.24 1.62 -81.89 -74.69 -73.08 -71.97 -51.94 ▁▇▁▁▁
zona 0 1.00 2.11 0.65 1.00 2.00 2.00 2.00 4.00 ▂▇▁▃▁
PESO_corr 0 1.00 755.69 1304.56 0.00 41.00 110.00 1079.66 14092.00 ▇▁▁▁▁
HORAS_REPOSO 2988 0.37 1629.72 832.08 12.00 1200.00 1200.00 2400.00 7200.00 ▇▅▁▁▁
N_PANOS 1833 0.61 226.33 280.51 1.00 1.00 2.00 420.00 900.00 ▇▁▃▁▁
order 0 1.00 2373.00 1369.91 1.00 1187.00 2373.00 3559.00 4745.00 ▇▇▇▇▇
CPUE 0 1.00 54.95 92.65 0.00 2.44 6.49 84.01 1067.68 ▇▁▁▁▁
CPUE2 2988 0.37 1.55 5.23 0.00 0.51 0.99 1.79 158.44 ▇▁▁▁▁

Datos 0

Are there missing values? y sirve para identificar los factores que se utilzaran en la estandarización.

dim(bitb3)
   [1] 4745   17
colSums(is.na(bitb3))
       N_TRIPULANTES POTENCIA_MOTOR_EQ               año               mes 
                   0               167                 0                 0 
                trim               dfp          lat_ctgr          lon_ctgr 
                   0                 0                51                51 
                zona         PESO_corr      TIPO_CARNADA      HORAS_REPOSO 
                   0                 0               122              2988 
             N_PANOS             order              CPUE             CPUE2 
                1833                 0                 0              2988 
        prof_med_cat 
                 939

Es probable que HORAS_REPOSOde la red no podamos usarlo como factor ni como variable de esfuerzo por lo exiguo del dato.

Variabilidad

asd1<-ggplot(data = bitb3 %>% 
               filter(zona!=4), 
             aes(y=dfp, x=zona,
                 group=zona)) +
  geom_boxplot()+
  labs(title = "Variabilidad en DFP",
       y = "DFP", x = "ZONA") + 
  facet_wrap(~ año, scales = "free")+
   theme_few()
asd1

asd2<-ggplot(data =bitb3%>%
               drop_na(prof_med_cat) %>% 
               filter(zona!=4), 
             aes(x=prof_med_cat,
                 y=lat_ctgr)) +
  geom_boxplot()+
  labs(title = "Variabilidad en distribución latitudinal",
       y = "Latitud", x = "Prof") + 
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  facet_wrap(~ año)
asd2

Normality and Homogenity factors

# Frequency polygon plot for catch
ggplot(bitb3 %>% 
         filter(zona!=4),
         aes(CPUE)) +
  geom_freqpoly(bins = 8) +
  labs(x = "Cpue", y = "Frequency") +
  facet_wrap(.~zona)+
  theme_few()+
  theme(panel.border = 
          element_rect(colour = "black", 
                       fill=NA, size = 1))

Shapiro-Wilk test for deviation from normality


shapiro.test(bitb3$CPUE)
   
    Shapiro-Wilk normality test
   
   data:  bitb3$CPUE
   W = 0.63039, p-value < 2.2e-16

Departure from normality…but not critical

Colinealidad entre factores

# COLLINEARITY
Coll <- c(
          "año",  "dfp" ,"lat_ctgr", "zona" ,
          "N_PANOS" ,    "CPUE", "prof_med_cat")

# Obtain summary using the ggpairs command from the GGally library
ggpairs(bitb3[,Coll], ggplot2::aes(alpha = 0.9, colour=prof_med_cat),
        lower = list(combo = "count"))+
  scale_fill_manual(values=c("red", "blue", "green",  "black")) +
   scale_colour_manual(values=c("red", "blue", "green", "black")) +
  theme_few()

# Nothing serious
#Calculate Variance Inflation Factor (VIF)
vifmodel <- round(vif(lm(PESO_corr ~ zona + dfp + año + mes,
                     data = bitb3)),2)
vifmodel
   zona  dfp  año  mes 
   1.34 1.29 1.06 1.01
barplot(vifmodel, main = "VIF Values", horiz = TRUE, col = "steelblue")

La multicolinealidad en el análisis de regresión se produce cuando dos o más variables predictoras están altamente correlacionadas entre sí, de modo que no proporcionan información única o independiente en el modelo de regresión.

Si el grado de correlación es lo suficientemente alto entre las variables, puede causar problemas al ajustar e interpretar el modelo de regresión.

La forma más común de detectar la multicolinealidad es mediante el factor de inflación de varianza (VIF), que mide la correlación y la fuerza de la correlación entre las variables predictoras en un modelo de regresión.

El valor de VIF comienza en 1 y no tiene límite superior. Una regla general para interpretar VIF es la siguiente:

Un valor de 1 indica que no hay correlación entre una variable predictora dada y cualquier otra variable predictora en el modelo. Un valor entre 1 y 5 indica una correlación moderada entre una variable predictora dada y otras variables predictoras en el modelo, pero esto a menudo no es lo suficientemente grave como para requerir atención. Un valor superior a 5 indica una correlación potencialmente grave entre una variable predictora dada y otras variables predictoras en el modelo. En este caso, es probable que las estimaciones de coeficientes y los valores p en el resultado de la regresión no sean fiables.

Interaccion

# Patterns in the variance? (Evidence for lack of homogeneity)
ggplot(bitb3, 
       aes(x = dfp, 
                  y = (CPUE))) +
  geom_jitter(shape = 16, 
              size = 2.5, 
              alpha = 0.3, 
              height = 0.25, 
              width = 0.5) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  facet_grid(bitb3$N_TRIPULANTES~bitb3$zona)+
  theme(panel.background = element_blank()) +
  theme(panel.border = element_rect(fill = NA, 
                                    size = 1)) +
  theme(strip.background = 
          element_rect(fill = "white",
                       color = "white",
                       size = 1)) +
  theme(text = element_text(size=13)) +
  theme_few()+
  xlab("Effort") + 
  ylab("CPUE") +
  ggtitle("variabless por Numero de Tripulantes y por zona")

Year x mes

ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~trim)

# No
# Year x habitat
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(zona~prof_med_cat)

# Perhaps
# Zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~trim)

  facet_grid(~zona)
   <ggproto object: Class FacetGrid, Facet, gg>
       compute_layout: function
       draw_back: function
       draw_front: function
       draw_labels: function
       draw_panels: function
       finish_data: function
       init_scales: function
       map_data: function
       params: list
       setup_data: function
       setup_params: function
       shrink: TRUE
       train_scales: function
       vars: function
       super:  <ggproto object: Class FacetGrid, Facet, gg>
# CPUE slope varies between habitats - interaction
# Tipo carnada zona
ggplot(bitb3,  aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few()+
  xlab("Points") + ylab("Catch") +
  facet_grid(TIPO_CARNADA~zona)

# Perhaps

Genero categorias y Factorizo

bitb3 <- bitb3 %>%
  mutate(POTENCIA_CAT = cut(POTENCIA_MOTOR_EQ, 
                            breaks = c(0, 200, 300, 400, 1000),
                      labels = c("Bajo", "Medio-Bajo", "Medio-Alto", "Alto")),
         LAT_CAT = cut(lat_ctgr, 
                              breaks = c(-20, -35, -40, -47),
                      labels = c("Norte", "Centro", "Sur")),
         N_TRIPULANTES=factor(N_TRIPULANTES),
         año=factor(año),
         zona=factor(zona),
         trim=factor(trim)) %>% 
  drop_na(zona,
          LAT_CAT,
          prof_med_cat,
          POTENCIA_CAT) %>%
  filter(CPUE>0) %>% 
  mutate(logCPUE=log(CPUE))

glimpse(bitb3)
   Rows: 3,347
   Columns: 20
   $ N_TRIPULANTES     <fct> 4, 8, 8, 8, 6, 7, 6, 7, 8, 8, 6, 8, 8, 9, 9, 9, 9, 9…
   $ POTENCIA_MOTOR_EQ <dbl> 250, 380, 320, 450, 380, 420, 370, 240, 320, 360, 45…
   $ año               <fct> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
   $ mes               <dbl> 7, 7, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11,…
   $ trim              <fct> 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
   $ dfp               <dbl> 13.833333, 11.741667, 15.086806, 13.729167, 18.31736
   $ lat_ctgr          <dbl> -40.66667, -39.96667, -40.50000, -33.91667, -40.6666
   $ lon_ctgr          <dbl> -74.50000, -74.33333, -74.50000, -72.58333, -74.6666
   $ zona              <fct> 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 3, 3, 3, 2, 2, 2, 2, 2…
   $ PESO_corr         <dbl> 1095.0568, 1256.6812, 1422.8584, 975.6000, 2026.6464
   $ TIPO_CARNADA      <chr> "2", "3", "3", NA, "1", "2", "2", "2", "3", "3", "3"
   $ HORAS_REPOSO      <dbl> 1400, 1800, 1200, 4800, 1200, 1200, 2400, 1200, 1200…
   $ N_PANOS           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
   $ order             <int> 8, 10, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 2…
   $ CPUE              <dbl> 79.160730, 107.027493, 94.311443, 71.060393, 110.640
   $ CPUE2             <dbl> 0.7821834, 0.6981562, 1.1857153, 0.2032500, 1.688872
   $ prof_med_cat      <fct> "(1e+03,1.5e+03]", "(500,1e+03]", "(1e+03,1.5e+03]",…
   $ POTENCIA_CAT      <fct> Medio-Bajo, Medio-Alto, Medio-Alto, Alto, Medio-Alto…
   $ LAT_CAT           <fct> Norte, Centro, Norte, Sur, Norte, Norte, Norte, Nort…
   $ logCPUE           <dbl> 4.371480, 4.673086, 4.546603, 4.263530, 4.706288, 4.

The data exploration showed:

  1. One significant outlier in catch
  2. Deviation from normality in response variable
  3. Possible departure from homogeneity
  4. Few zeros in the response variable
  5. No collinearity
  6. No imbalance
  7. Possible season x habitat interaction

Chequeo distribuciones


bcpue <- ggplot(bitb3,
                aes(CPUE))+
  coord_flip()+
  geom_boxplot()+
  theme_few()


hcpue <- ggplot(bitb3,
                aes(CPUE))+
  geom_histogram(bins=15, fill=2)+
  theme_few()


ggarrange(bcpue, hcpue, ncol = 2)  


blcpue <- ggplot(bitb3,
                aes(logCPUE))+
  coord_flip()+
  geom_boxplot()+
  theme_few()


hlcpue <- ggplot(bitb3,
                aes(logCPUE))+
  geom_histogram(bins=15, fill=2)+
  theme_few()


ggarrange(blcpue, hlcpue, ncol = 2)  

#Chequeo visual de normalidad y dispersion
#-------------------------------------------------------------------------------
#CPUE y LOGCPUE
par(mfrow=c(2,2))

qqnorm(bitb3$CPUE)   
qqline(bitb3$CPUE,col='red') #
boxplot(CPUE ~ año,data=bitb3,na.rm=T,#
        main="Boxplot nominal CPUE by year", ylab="CPUE",col='light grey',boxwex=0.65)

qqnorm(bitb3$logCPUE)   
qqline(bitb3$logCPUE,col='red') #
boxplot(logCPUE ~ año,data=bitb3,na.rm=T,#
        main="Boxplot nominal LogCPUE by year", ylab="LogCPUE",col='light grey',boxwex=0.65)

Modelos GLM

#Predictores (Factores Principales)
gauss01 = glm(formula= logCPUE ~  año,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss02 = glm(formula= logCPUE ~  año+trim,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss03 = glm(formula= logCPUE ~  año+trim+prof_med_cat,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss04 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss05 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss06 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+TIPO_CARNADA,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss07 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
gauss08 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona:prof_med_cat,
              family = gaussian(link="identity"), 
              data = bitb3,na.action=na.exclude)
rmodelo01 = c(AIC(gauss01),(gauss01$null.deviance-gauss01$deviance)/gauss01$null.deviance)
rmodelo02 = c(AIC(gauss02),(gauss02$null.deviance-gauss02$deviance)/gauss02$null.deviance)
rmodelo03 = c(AIC(gauss03),(gauss03$null.deviance-gauss03$deviance)/gauss03$null.deviance)
rmodelo04 = c(AIC(gauss04),(gauss04$null.deviance-gauss04$deviance)/gauss04$null.deviance)
rmodelo05 = c(AIC(gauss05),(gauss05$null.deviance-gauss05$deviance)/gauss05$null.deviance)
rmodelo06 = c(AIC(gauss06),(gauss06$null.deviance-gauss06$deviance)/gauss06$null.deviance)
rmodelo07 = c(AIC(gauss07),(gauss07$null.deviance-gauss07$deviance)/gauss07$null.deviance)
rmodelo08 = c(AIC(gauss08),(gauss08$null.deviance-gauss08$deviance)/gauss08$null.deviance)
resultados = rbind(rmodelo01,rmodelo02,rmodelo03,rmodelo04,rmodelo05,rmodelo06,
                   rmodelo07,rmodelo08);resultados
check_model(gauss07)

Table with significance level (***).


tab_model(gauss01, 
          gauss02,
          gauss03, 
          gauss04,
          gauss05,
          gauss06,
          gauss07,
          gauss08,
          p.style = "stars")
model_performance(gauss07)
   # Indices of model performance
   
   AIC       |      AICc |       BIC |    R2 |  RMSE | Sigma
   ---------------------------------------------------------
   12846.652 | 12847.141 | 13017.894 | 0.335 | 1.635 | 1.642

Table comparing performance model.

compare_performance(gauss01, 
          gauss02,
          gauss03, 
          gauss04,
          gauss05,
          gauss06,
          gauss07,
          gauss08,
          rank = TRUE, 
          verbose = FALSE)
   # Comparison of Model Performance Indices
   
   Name    | Model |    R2 |  RMSE | Sigma | Performance-Score
   -----------------------------------------------------------
   gauss08 |   glm | 0.341 | 1.627 | 1.635 |           100.00%
   gauss07 |   glm | 0.335 | 1.635 | 1.642 |            97.35%
   gauss06 |   glm | 0.334 | 1.643 | 1.650 |            95.42%
   gauss05 |   glm | 0.299 | 1.679 | 1.685 |            82.38%
   gauss04 |   glm | 0.264 | 1.720 | 1.726 |            68.13%
   gauss03 |   glm | 0.217 | 1.774 | 1.778 |            49.32%
   gauss02 |   glm | 0.102 | 1.900 | 1.904 |             4.17%
   gauss01 |   glm | 0.091 | 1.912 | 1.915 |             0.00%

Regarding this result, we can select best model and check some assumtion like:

  • Colinealidad
  • Data points influyentes
  • Homoscedasticidad
  • Normalidad de los residuos
  • Independencia

COMPOSICIONES DE TALLAS

La información contenida en las estructiras de tallas son consideradas una de las piezas mas importantes en la ciencia pesquera (Canales, Punt, & Mardones, 2021; Hordyk, Ono, Sainsbury, Loneragan, & Prince, 2014; Rudd & Thorson, 2018). Sin embargo, esta debe ser validada y confirmada para no violar supuestos referidos a la repreentatividad del dato y su relación con la dinámica poblacional.

Ahora nos disponemos a explorar esta fuente de información.

Primero identificamos la estructura de la base;

glimpse(long)
   Rows: 15,831
   Columns: 24
   $ ID                  <chr> "COD_BARCO=730881 AND TO_CHAR(FECHA_HORA_RECALADA,…
   $ COD_BARCO           <dbl> 730881, 730881, 730881, 730881, 730881, 730881, 73…
   $ FECHA_HORA_RECALADA <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
   $ FECHA_HORA_ZARPE    <dttm> 2021-05-10 11:00:00, 2021-05-10 11:00:00, 2021-05…
   $ COD_PESQUERIA       <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50…
   $ PUERTO_RECALADA     <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 3,…
   $ NRO_FORMULARIO      <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,…
   $ REGION              <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1,…
   $ FECHA_LANCE         <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
   $ NUMERO_LANCE_EX     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
   $ LATITUD             <dbl> 243900, 243900, 243900, 243900, 243900, 243900, 24…
   $ LONGITUD            <dbl> 714400, 714400, 714400, 714400, 714400, 714400, 71…
   $ NUMERO_EJEMPLARES   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ PESO_TOTAL_MUESTRA  <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ NUMERO_CAJA         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
   $ COD_ESPECIE         <dbl> 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37…
   $ ORIGEN_MUESTRA      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3,…
   $ FECHA_MUESTREO      <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ N_TOTAL_INDIV       <dbl> 85, 85, 85, 85, 85, 85, 85, 84, 84, 84, 84, 84, 23…
   $ LONGITUD_DESCARTE   <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ LONGITUD_MUESTRA    <dbl> 115, 125, 129, 130, 143, 152, 157, 100, 108, 110, …
   $ SEXO                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
   $ N_INDIVIDUOS        <dbl> 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 1, 1,…
   $ año                 <dbl> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 20…
colSums(is.na(long))
                    ID           COD_BARCO FECHA_HORA_RECALADA    FECHA_HORA_ZARPE 
                  2997                   0                   0                   0 
         COD_PESQUERIA     PUERTO_RECALADA      NRO_FORMULARIO              REGION 
                     0                   0                   0                   0 
           FECHA_LANCE     NUMERO_LANCE_EX             LATITUD            LONGITUD 
                     0                   0                4425                4425 
     NUMERO_EJEMPLARES  PESO_TOTAL_MUESTRA         NUMERO_CAJA         COD_ESPECIE 
                 15733               15831                   0                   0 
        ORIGEN_MUESTRA      FECHA_MUESTREO       N_TOTAL_INDIV   LONGITUD_DESCARTE 
                     0               15831                   0               15831 
      LONGITUD_MUESTRA                SEXO        N_INDIVIDUOS                 año 
                     0                   0                   2                   0

Primero selecciono las columnas de interes y luego genero la expansión de LONGITUD_MUESTRA a N_INDIVIDUOS. Expand frecuency data related length, in this case N_INDIVIDUOS column have frecuency that we need expand to whole data frame.

Selecciono variables de interes

long1 <- long %>% 
  dplyr::select(6,8,11,12,21,22,23,24)
dim(long1)
   [1] 15831     8
long2 <- long %>% 
  drop_na(N_INDIVIDUOS) %>% 
  type.convert(as.is = TRUE) %>% 
  uncount(N_INDIVIDUOS)

dim(long2)
   [1] 34299    23

legend.labels <- c('Indet.', 'Macho' , 'Hembra')
nbco <- ggplot(long2, aes(x=LONGITUD_MUESTRA, y = as.factor(año),
                         fill=as.factor(SEXO)))+
  geom_density_ridges(stat = "density_ridges", bins = 20, 
                      scale = 2, draw_baseline = FALSE,
                      alpha=0.5)+
  facet_wrap(.~REGION, ncol=5) +
  geom_vline(xintercept = 110, color = "red")+
  scale_fill_viridis_d(option = "C",
                       name="SEXO",
                       labels=legend.labels)+
  scale_y_discrete(breaks = seq(from = 2004, to = 2022, by = 2))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  theme_minimal()+
  xlab("Longitud (cm.)")+
  ylab("")
#scale_x_discrete((limits = rev(levels(talla2021$ANO_ARR))))+
nbco

Prepara los vectores para sumar a los .dat del modelo si los necesito.

long2$LONG_CAT <- as.numeric(as.character(cut(x = long2$LONGITUD_MUESTRA, 
                                              breaks = seq(10,220,2), 
                                              labels = seq(10,218,2),
                                              right = FALSE)))
LONGT <- table(long2$año, long2$LONG_CAT)
LONGT
         
           10  12  38  40  42  44  48  50  52  54  56  58  60  62  64  66  68  70
     2004   2   1   0   0   0   0   0   0   0   0   0   3  14  17   3  13  14  23
     2005   0   0   0   0   0   0   0   0   0   0   1   0   1   2   5  10  11  17
     2006   0   0   0   0   0   0   0   0   0   2   2   3   2   4   6   6   8  12
     2007   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     2008   0   0   0   0   0   0   0   0   1   0   0   0   2   3   2   3  10   3
     2009   0   0   0   0   0   0   0   0   0   0   1   0   5   2   4   8  25  23
     2010   0   0   0   0   0   0   0   0   0   1   1   2   1   5   7  11  11  16
     2011   0   0   0   1   2   1   1   1   0   1   2   2   4   9  12   8  17  21
     2012   0   0   0   0   0   0   0   0   0   0   0   0   1   5   2  20  20  44
     2013   0   0   0   0   0   0   0   0   0   1   0   0   0   1   0   0  12  10
     2014   0   0   1   0   0   0   1   0   0   0   1   2  27  21  41  85 112 120
     2015   0   0   0   0   0   0   0   0   0   0   0   1   6  16  20  41  65  87
     2016   0   0   0   0   0   0   0   0   0   0   1   4  10  17  47  56  90 140
     2017   0   0   1   0   0   0   0   0   0   3   9  19  60  59 122 117 111 178
     2018   0   0   0   0   0   0   0   0   0   6   8  25  27  35  46  61  47  68
     2019   0   0   0   0   0   0   0   0   0   1   2   4  18  42  81 118 139 202
     2020   0   0   0   0   0   0   0   0   0   0   0   0   4   5  11  17  11  20
     2021   0   0   0   0   0   0   0   0   1   0   0   0   3   4  10  10  18  41
     2022   0   0   0   0   0   0   0   0   2   1   5  13  24  32  48  83  88 139
         
           72  74  76  78  80  82  84  86  88  90  92  94  96  98 100 102 104 106
     2004  22  22  40  28  33  50  31  62  54  72  47  37  41  43  52  21  17  22
     2005  19  19  26  21  21  31  31  22  31  25  22  18  15  22  10  17   9   8
     2006  11  13  17  20  26  13  12  11   9   9   7   3  12   6   6   9   8   7
     2007   0   1   0   0   1   1   1   2   0   2   0   0   2   0   0   0   0   0
     2008  17   9  17  22  12  15  13  17  16  10  12  10   9  10   6   7   6  10
     2009  39  52  83  82  84 113 105  94  97  82  98  71  54  50  38  28  43  15
     2010  14  22  31  36  21  34  36  41  31  57  36  28  28  20  30  20  11  15
     2011  11  20  25  21  11  18  13  12  18   9  12  10  13  13   8   9   5   8
     2012  63  67  60  39  61  32  29  24  18  27  34  14  12   6  13   9  14   5
     2013  29  29  37  61  99  82  95 112  62 115  55  70  68  41  44  42  34  43
     2014 186 188 216 230 240 252 243 241 224 190 166 124 103 101  68  70  51  36
     2015 116 165 186 186 235 248 253 268 195 240 237 193 179 142 139  86  81  47
     2016 147 192 176 223 278 280 270 270 228 239 234 228 166 154 154 140  99  88
     2017 159 171 146 152 158 160 196 178 147 190 162 184 180 160 217 131 139  94
     2018  46  50  33  28  35  25  23  43  33  43  32  45  44  41  37  36  30  41
     2019 217 190 165 124 102  80  54  31  28  39  24  33  26   8  43  23  26  27
     2020  29  30  28  19  21  23  13   1   6   8   3   4   1   7   4   4   2   7
     2021  54  77 104 115 179 143 159 128 107 114  75  61  48  30  36  25  19  14
     2022 173 203 224 270 323 372 383 376 303 335 285 252 192 119 123 119  76  56
         
          108 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142
     2004  15  16  11   9  14  11  21   7   6   2   4   8   9   5   4   2   6   3
     2005  12  10   8   7   5   3   3   3   2   5   1   2   1   0   2   1   2   0
     2006   6   2   5   5   2   2   4   4   3   2   1   0   1   3   1   0   0   2
     2007   0   0   0   1   0   0   0   0   0   0   1   0   0   0   0   0   0   0
     2008   7   6   7  11   6   6   2   4   2   4   5   4   2   2   3   1   3   3
     2009  23  25  18  19  15  12  13   9  12   8   9  10   4   8  11  10  15  11
     2010  16  13   6   8   5   7   9   4   4   5   0   1   1   1   2   4   1   2
     2011   7  15  15   4   1   7   2   4   9   4   4   2   1   1   2   3   2   1
     2012   8  10   3   3   2   3   2   0   1   1   0   0   1   0   1   0   2   0
     2013  31  35  21  21  15  12  21   8   7   7   6   9   5   4   3   4   5   0
     2014  46  35  38  38  41  35  41  29  26  13   9  17  12   8   9   7   9  11
     2015  52  52  39  34  44  41  33  19  27  29  16  18  18   8  15  12  15   5
     2016  45  58  42  30  30  29  28  22  20  23  11  18  13   9   5   6   7   5
     2017  77 139  60  78  45  37  55  19  24  20  11  26  18  22  17   5  12  10
     2018  30  25  26  22  14  10  21  14   4   9   6  10   2   9   9   5   5   7
     2019  24  26  20  22   7  17  22  15  14  11   7   5   1   7   8   1   6   3
     2020   3   3   2   5   4   4   4   0   2   0   4   2   1   0   1   1   1   0
     2021  13  26  34  20  29  21  32  21  27  27  20  17  21  12  20  15  17  13
     2022  33  52  22  21  35  29  14  22  26  21  37  24  21  27  25  14  19  20
         
          144 146 148 150 152 154 156 158 160 162 164 166 168 170 172 174 176 178
     2004   0   3   2   1   0   0   0   0   0   0   0   1   0   0   0   0   0   0
     2005   1   1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     2006   0   0   1   0   0   1   0   0   0   0   1   0   0   0   0   0   0   0
     2007   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     2008   1   3   5   1   3   1   3   3   0   3   0   0   0   0   0   0   0   0
     2009   4   7   9  13  10   2   8  10   5   5   2   5   5   0   1   1   2   0
     2010   5   0   2   3   0   2   1   2   4   0   0   0   1   1   1   0   0   0
     2011   1   0   2   1   0   0   0   0   1   0   0   1   0   1   0   0   0   0
     2012   1   0   0   1   0   1   0   0   0   0   0   1   0   0   2   0   0   1
     2013   5   4   2   0   0   2   1   0   0   0   2   1   0   0   0   1   0   0
     2014   5   5   6   1   1   2   1   1   3   2   2   3   1   1   0   0   1   0
     2015  10   6  10   9   5   5   3   4   2   2   0   3   3   3   2   1   1   0
     2016   7   7   4   4   4   5   0   1   0   4   3   5   1   2   0   0   3   0
     2017   8   9   7   9   8   4   0   5   6   6   4   5   0   3   2   1   2   0
     2018   5   1   4   6   2   1   1   5   3   1   1   0   1   2   1   3   0   0
     2019   3   4   2   4   1   1   3   0   2   0   2   0   0   0   0   0   0   0
     2020   2   0   0   0   1   0   0   0   2   0   0   0   0   2   0   0   0   0
     2021  11   7   5  13   7   2   5   5   8   6   2   2   4   5   0   0   0   1
     2022  19  18  16  22  17   5  10   7  13   9   6   5   7   5   3   4   3   2
         
          180 182 184 186 188 190 192 200 202
     2004   1   0   0   0   0   0   0   0   0
     2005   0   0   0   0   0   0   0   0   0
     2006   0   0   0   0   0   0   0   0   0
     2007   0   0   0   0   0   0   0   0   0
     2008   0   0   0   0   0   0   0   0   0
     2009   0   0   1   0   0   0   0   0   0
     2010   0   1   0   0   0   0   0   0   0
     2011   0   0   0   0   0   0   0   0   0
     2012   0   0   0   0   0   0   0   0   0
     2013   0   0   0   0   0   0   0   0   0
     2014   0   0   0   0   0   0   1   0   0
     2015   0   0   0   0   0   0   0   0   0
     2016   0   0   1   0   1   0   0   0   0
     2017   0   1   0   0   0   0   0   1   0
     2018   0   0   0   0   0   0   0   0   1
     2019   0   0   0   0   0   0   0   0   0
     2020   0   0   0   0   0   0   0   0   0
     2021   0   0   0   0   0   0   0   0   0
     2022   0   0   3   1   0   1   0   0   0
write.csv()

BITACORA FIPA 96- 14


bitfipa<- bit8696 %>% 
  mutate(CPUE = CAPTURA/DFP) %>% 
  drop_na(CPUE)


cpuefipa <- ggplot(bitfipa %>%
                     filter(ANO<1997)  %>%
                  group_by(ANO) %>% 
                  summarise(CPUEF=mean(CPUE)),
                aes(ANO, CPUEF))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1986, to = 2022, by = 1))+
  geom_smooth(method = 'lm', 
              colour = 'blue', 
              size = 1.5)+
  # stat_regline_equation(label.x=1993, label.y=500)+
  # stat_cor(aes(label=..rr.label..), label.x=1993, label.y=450)+
  #facet_wrap(~AREA)+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="CPUE (kg/dfp)",
       x="")
cpuefipa

Guardo datos FIPA 96

cpuenomfipa <- bitfipa %>%
  filter(ANO<1997)  %>%
  group_by(ANO) %>% 
  summarise(CPUEF=mean(CPUE))
write_csv(cpuenomfipa, "CPUE_FIPA.csv")

MAPAS

Lo primero es transformar los datos de lon_ctgr y lat_ctgr a formarto geom. A su vez, saco las NA de las coord

transformar los datos en un sf object

codmap <- st_as_sf(bitb2 %>% 
                    drop_na(lon_ctgr) %>% 
                    drop_na(lat_ctgr), 
                   coords = c("lon_ctgr", "lat_ctgr"),  crs = 4326)

Ahora genero los mapas de Chile con un raster.

chile <- raster::getData("GADM", country = "CHL", level = 0)
chile1<-fortify(chile)
chilemap <- ggplot()+
   geom_polygon(data=chile, aes(x=long, y=lat, group=group),
              fill="lightblue",color="grey20", size=0.15)+
   coord_sf(crs = st_crs(4326),
            xlim = c(-80, -65),
            ylim = c(-58, -15)) +
          theme_bw()
 chilemap

# Aca veo los nombres
chile@data$NAME_1
   NULL

Luego genero los bordes sobre los cuales haré la grilla.

e <- extent(-80,-65,-58,-15)
rc <- crop(chile, e)

# la proyección adecuada en lat long
tcrs <- CRS("+init=epsg:4326")
transformed_points <- spTransform(rc, tcrs)
# para dejarlo en formato geom_sf
chile2 <- st_as_sf(transformed_points) 

Se estructura la grilla en el objeto raster de chile1

grid <- chile2 %>%
  st_make_grid(cellsize = c(1,1)) %>% # para que quede cuadrada
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>% # objeto en spatial feature
  mutate(cellid = row_number())

Pongo los datos codmap en la grilla. Aca solo elegí dfpy CPUE pero se ùeden resumir otros

joindat <- grid %>%
  st_join(codmap) %>% 
  group_by(cellid, año) %>% 
  summarise(CPUEM = mean(CPUE),
            DFPM =mean(dfp)) # por ejemplo, plotear la captura "CAPTURA_1"

Mapa Esfuerzo

here we can viz diferent variables

## Plot final
mas1 <- ggplot() +
  geom_sf(data=joindat %>% 
             filter(!is.na(DFPM)), aes(fill = DFPM),
          color=NA) +
  scale_fill_viridis_b(option="E",
                       direction=-1, name="DFP")+
  geom_sf(data = grid,  fill=NA, color=NA) +
  geom_sf(data = chile2, color="grey", fill="white") +
  coord_sf() +
  geom_hline(yintercept = -47, color = "red")+
  scale_alpha(guide="none")+
  facet_wrap(~año, ncol=6)+
  scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides(colour = guide_legend()) +
  theme_bw()+
  theme(panel.background = element_rect(fill = 'white'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())
mas1

Mapa CPUE

## Plot final
mas2 <- ggplot() +
  geom_sf(data=joindat %>% 
             filter(!is.na(CPUEM)), aes(fill = CPUEM),
          color=NA) +
  scale_fill_viridis_b(option="G",
                       direction=-1, name="CPUE")+
  geom_sf(data = grid,  fill=NA, color=NA) +
  geom_sf(data = chile2, color="grey", fill="white") +
  coord_sf() +
  geom_hline(yintercept = -47, color = "red")+
  scale_alpha(guide="none")+
  facet_wrap(~año, ncol=6)+
  scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides(colour = guide_legend()) +
  theme_bw()+
  theme(panel.background = element_rect(fill = 'white'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())
mas2

Falta componer mapas de tallas medias y alguna otra variable de interés.

CONCLUSION

En función de los datos analizados, el modelo propuesto es un enfoque de producción global estado espacio descrito por Pedersen & Berg (2017) y Acom (2015).

REFERENCIAS

Acom, I. C. M. (2015). ICES WKLIFE V REPORT 2015 Methodologies based on Life-history Traits , Development of Quantitative Assessment Exploitation Characteristics and other Report of the Fifth Workshop on the Relevant Parameters for Data-limited Stocks ( WKLIFE V ) Lisbon , Port. (October), 5–9.
Canales, C. M., Punt, A. E., & Mardones, M. (2021). Can a length-based pseudo-cohort analysis (LBPA) using multiple catch length-frequencies provide insight into population status in data-poor situations? Fisheries Research, 234(October 2020), 105810. https://doi.org/10.1016/j.fishres.2020.105810
Hordyk, A., Ono, K., Sainsbury, K., Loneragan, N., & Prince, J. (2014). Some explorations of the life history ratios to describe length composition, spawning-per-recruit, and the spawning potential ratio. ICES Journal of Marine Science, 72(1), 204–216. https://doi.org/10.1093/icesjms/fst235
Pedersen, M. W., & Berg, C. W. (2017). A stochastic surplus production model in continuous time. Fish and Fisheries, 18(2), 226–243. https://doi.org/10.1111/faf.12174
Rudd, M. B., & Thorson, J. T. (2018). Accounting for variable recruitment and fishing mortality in length-based stock assessments for data-limited fisheries. Canadian Journal of Fisheries and Aquatic Sciences, 75(7), 1019–1035. https://doi.org/10.1139/cjfas-2017-0143